import numpy as np

def f1(a, r, d):
    return np.cos(a) - (1 - d**2 / (2 * r**2))

def f2(a, r, L):
    return L - a * r

def jacobian(a, r, d):
    return np.array([[-np.sin(a), d**2 / r**3], [r, a]])

# 给定参数
d = 140
L = 156

# 初始猜测
a, r = 1, 1

# 容差和最大迭代次数
tol = 1e-6
max_iter = 100

for _ in range(max_iter):
    F = np.array([f1(a, r, d), f2(a, r, L)])
    J = jacobian(a, r, d)
    delta = np.linalg.solve(J, -F)
    a, r = a + delta[0], r + delta[1]
    
    if np.linalg.norm(F) < tol:
        break

print(f"Solutions: a = {a}, r = {r}")